1 Introduction

This file was last update on 2020-09-03

2 Sampling effort

No. of visits to each transect/site.

abn%>%
  group_by(Study, Taxa, Transect)%>%
  summarise(n = length(unique(Date)))%>%
  spread(Transect, n)
Study Taxa T1A T1B T2 T3 T4 T5 T6
1998 - 2000 Birds 19 18 7 10 4 1 1
1998 - 2000 Butterflies 14 7 4 9 4 1 1
2016 - 2017 Birds 23 23 22 19 16 22 17
2016 - 2017 Butterflies 23 20 20 22 18 21 15

Uneveness in sampling effort across sites in studies may bias estimates of change.

3 Change in diversity across years

3.1 Species richness

3.1.1 Across the two studies

abn%>%
  group_by(Taxa, Study)%>%
  summarise(richness = length(unique(Specific.Epithet)),
            N = sum(Abundance, na.rm = T))
Taxa Study richness N
Birds 1998 - 2000 70 1007
Birds 2016 - 2017 105 2021
Butterflies 1998 - 2000 45 515
Butterflies 2016 - 2017 66 2014

More species of both taxa were encountered in current study compared to the previous study. Difference in sampling effort is apparent accross two studies.

3.1.2 Across sites in the two studies

No. of species encountered across sites sampled in the two studies.

abn%>%
  left_join(transects, by = c("Transect"))%>%
  group_by(Taxa, Study, Transect)%>%
  summarise(richness = length(unique(Specific.Epithet)))%>%
  spread(Transect, richness)
Taxa Study T1A T1B T2 T3 T4 T5 T6
Birds 1998 - 2000 47 37 28 36 29 4 11
Birds 2016 - 2017 66 59 53 47 31 38 28
Butterflies 1998 - 2000 27 22 8 37 24 5 3
Butterflies 2016 - 2017 41 32 36 44 31 39 27

Fewer species were encountered in sites with lower effort in the previous study, i.e., T5 and T6.

3.2 Shannon diversity

# Calculating diversity in each site over years

div_year <- abn%>%
  group_by(Study, Taxa, Transect,  Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  spread(Specific.Epithet, n, fill = c(0))%>%
  group_by(Study, Taxa, Transect)%>%
  nest()%>%
  mutate(H = map_dbl(data, ~vegetarian::d(.[-1], lev = "alpha", q = 1)))%>%
  dplyr::select(-data)

# Summary

div_year%>%
  group_by(Taxa, Study)%>%
  skimr::skim(H)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

Variable type: numeric

Taxa Study mean sd
Birds 1998 - 2000 15.39 7.79
Birds 2016 - 2017 24.88 5.74
Butterflies 1998 - 2000 11.06 6.57
Butterflies 2016 - 2017 20.00 3.01 
# t test

div_year%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(H ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()
Taxa 1998-2001 2016-2017 parameter statistic p.value
Birds 15.39287 24.87583 11.031229 -2.593073 0.0249503
Butterflies 11.05543 19.99721 8.410531 -3.273453 0.0105504 
# plot

div_year%>%
  ggplot(aes(Taxa, H, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Effective number of species")+
  scale_fill_brewer(palette = "Greys")

ggsave(last_plot(), filename = "./Figures and Tables/figure2.png", height = 6, width = 8, dpi = 300)

The first order diversity of both taxa has increased significantly across the two surveys.

4 Change in composition of species accross years

library(vegan)

# input for vegan

veg_birds <- abn%>%
  filter(Taxa == "Birds")%>%
  group_by(Study, Transect, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = paste(Study, "_", Transect))%>%
  group_by(sample)%>%
  #mutate(N = sum(n, na.rm = T))%>%
  #filter(N > 0)%>%#removing samples where no species were detected
  dplyr::select(sample, Specific.Epithet, n)%>%
  ungroup()%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

veg_butterfly <- abn%>%
  filter(Taxa == "Butterflies")%>%
  group_by(Study, Transect, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = paste(Study, "_", Transect))%>%
  group_by(sample)%>%
  #mutate(N = sum(n, na.rm = T))%>%
  #filter(N > 0)%>%#removing samples where no species were detected
  dplyr::select(sample, Specific.Epithet, n)%>%
  ungroup()%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

# sample info

samp_birds <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  filter(Taxa == "Birds")%>%
  dplyr::select(Study, Taxa, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))

samp_butterfly <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  filter(Taxa == "Butterflies")%>%
  dplyr::select(Study, Taxa, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))
# Similarity matrix by studies

sim_birds <- abn%>%
  filter(Taxa == "Birds")%>%
  group_by(Study, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = Study)%>%
  ungroup()%>%
  dplyr::select(sample, Specific.Epithet, n)%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

sim_butterfly <- abn%>%
  filter(Taxa == "Butterflies")%>%
  group_by(Study, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = Study)%>%
  ungroup()%>%
  dplyr::select(sample, Specific.Epithet, n)%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

# bray - curtis

data_frame(
  Taxa = c("Birds", "Butterflies"),
  Dissimilarity = c(as.numeric(vegdist(sim_birds, method = "bray")),
                    as.numeric(vegdist(sim_butterfly, method = "bray")))
)
Taxa Dissimilarity
Birds 0.5495376
Butterflies 0.6710162 
# ordination

NMDS_birds <- metaMDS(veg_birds, distance = "bray", k = 2, model = "local")

NMDS_butterfly <- metaMDS(veg_butterfly, distance = "bray", trymax = 1000, k = 2, model = "local")


# Extracting scores and adding habitat vars

ord_birds <- data.frame(scores(NMDS_birds))%>%
  rownames_to_column("sample")%>%
  left_join(samp_birds, "sample")

ord_butterfly <- data.frame(scores(NMDS_butterfly))%>%
  rownames_to_column("sample")%>%
  left_join(samp_butterfly, "sample")

ord <- bind_rows(ord_birds, ord_butterfly)
#Testing change in composition of birds

adonis(veg_birds ~ Study, data = samp_birds)
## 
## Call:
## adonis(formula = veg_birds ~ Study, data = samp_birds) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Study      1   0.77547 0.77547  3.9632 0.24827  0.001 ***
## Residuals 12   2.34798 0.19567         0.75173           
## Total     13   3.12345                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Testing change in composition of butterflies

adonis(veg_butterfly ~ Study, data = samp_butterfly)
## 
## Call:
## adonis(formula = veg_butterfly ~ Study, data = samp_butterfly) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Study      1    0.9191 0.91908   4.054 0.25253  0.001 ***
## Residuals 12    2.7205 0.22671         0.74747           
## Total     13    3.6396                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# plotting

ord%>%
  ggplot(aes(NMDS1, NMDS2, col = Study, fill = Study))+
  stat_ellipse(level = 0.95, geom = "polygon", linetype = "dashed", size = 1, alpha = 0.3)+
  geom_point(aes(shape = Habitat), size = 3)+
  guides(shape = guide_legend(nrow = 2, byrow = T))+
  scale_color_brewer(palette = "Set1")+
  scale_fill_brewer(palette = "Set1")+
  facet_wrap(~Taxa)+
  theme(legend.box = "vertical")

ggsave(last_plot(), filename = "./Figures and Tables/figure3.png", height = 6, width = 8, dpi = 300)

Both bird and butterfly communities compositions are significantly different compared to the previous study.

5 Change in niche width of community

5.1 Trophic Guilds

# compiling trophic guild data

host <- read.csv("./Data/Butterflies_host plant_280820.csv") #butterfly host plant data

diet <- read.csv("./Data/Birds_diet.csv") #bird diet data

trophic <- host%>%
  full_join(diet)%>%
  dplyr::select(Specific.Epithet, Guild)%>%
  distinct()

# Calculating relative proportion of guilds

guilds <- abn%>%
  left_join(trophic, by = c("Specific.Epithet"))%>%
  group_by(Study, Taxa, Transect,  Guild)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  #group_by(Study, Taxa, Transect)%>%
  mutate(N = sum(n))%>%
  #filter(Guild == "Generalist" | Guild == "Omnivore")%>%
  mutate(prop.gen = n/N)%>%
  filter(Guild != "",
         !is.na(Guild))
# summary

guilds%>%
  group_by(Taxa, Guild, Study)%>%
  skimr::skim(prop.gen)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Guild, Study, mean, sd)

Variable type: numeric

Taxa Guild Study mean sd
Birds Carnivore 1998 - 2000 0.15 0.09
Birds Carnivore 2016 - 2017 0.04 0.03
Birds Frugivore 1998 - 2000 0.25 0.08
Birds Frugivore 2016 - 2017 0.28 0.05
Birds Grainivore 1998 - 2000 0.09 0.04
Birds Grainivore 2016 - 2017 0.04 0.02
Birds Insectivore 1998 - 2000 0.17 0.03
Birds Insectivore 2016 - 2017 0.36 0.07
Birds Omnivore 1998 - 2000 0.36 0.12
Birds Omnivore 2016 - 2017 0.29 0.06
Butterflies Generalist 1998 - 2000 0.22 0.09
Butterflies Generalist 2016 - 2017 0.18 0.07
Butterflies Grass Specialist 1998 - 2000 0.11 0.09
Butterflies Grass Specialist 2016 - 2017 0.17 0.08
Butterflies Herb Specialist 1998 - 2000 0.11 0.04
Butterflies Herb Specialist 2016 - 2017 0.21 0.12
Butterflies Liana Specialist 1998 - 2000 0.01 0.00
Butterflies Shrub Specialist 1998 - 2000 0.32 0.19
Butterflies Shrub Specialist 2016 - 2017 0.14 0.04
Butterflies Tree Specialist 1998 - 2000 0.38 0.09
Butterflies Tree Specialist 2016 - 2017 0.30 0.04 
# t - test

guilds%>%
  filter(Guild != "Liana Specialist")%>% #removing due to lack of obs
  complete(nesting(Study, Taxa), Transect, Guild, fill = list(prop.gen = 0))%>%
  group_by(Taxa, Guild)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(prop.gen ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()
Taxa Guild 1998-2001 2016-2017 parameter statistic p.value
Birds Carnivore 0.1507973 0.0263170 63.83105 10.037830 0.0000000
Birds Frugivore 0.2455696 0.2786997 82.30864 -2.535466 0.0131216
Birds Grainivore 0.0943404 0.0396365 79.81348 9.738085 0.0000000
Birds Insectivore 0.1430346 0.3600924 87.58662 -15.598448 0.0000000
Birds Omnivore 0.3647138 0.2936409 71.97471 3.934022 0.0001905
Butterflies Generalist 0.1590606 0.1817130 46.39068 -1.011897 0.3168362
Butterflies Grass Specialist 0.0792936 0.1697149 68.44700 -4.838736 0.0000078
Butterflies Herb Specialist 0.0647845 0.2118724 74.95577 -7.332093 0.0000000
Butterflies Shrub Specialist 0.3178457 0.1356982 52.10668 7.205350 0.0000000
Butterflies Tree Specialist 0.3761303 0.3010015 63.30807 5.941932 0.0000001 
# plot
guilds%>%
  ggplot(aes(reorder(Guild, prop.gen), prop.gen, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5, position = position_dodge(preserve = "single"))+
  labs(y = "Relative Proportion", x = "Trophic Guild")+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  scale_fill_brewer(palette = "Greys")+
  facet_wrap(~Taxa, scales = "free_x", ncol = 1)

ggsave(last_plot(), filename = "./Figures and Tables/figure4.png", height = 9, width = 8, dpi = 300)

Carnivorous, Grainivorous birds have reduced in the comunity, compared to insectivorous birds which have increased. Similarly, Shrub specialist and tree specialist butterflies have reduced where grass specialists and herb specialist have increased. Suprsingly, generalist species of both taxa appear to be unchanged.

6 Change in niche specialisation

-Can we incorporate abunndance into community specialisation index?

6.1 Trophic niche

# trophic specialisation based on Juliard et al. 2006

butterflies_TSI <- host%>%
  mutate(H = length(unique(Hostplant.Family)))%>%
  group_by(Specific.Epithet)%>%
  mutate(h = length(unique(Hostplant.Family)))%>%
  ungroup()%>%
  mutate(TSI = sqrt((H/h)-1))

birds_TSI <- diet%>%
  mutate(H = length(unique(Prey)))%>%
  group_by(Specific.Epithet)%>%
  mutate(h = length(unique(Prey)))%>%
  ungroup()%>%
  mutate(TSI = sqrt((H/h)-1))

CSI.T <- abn%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
  group_by(Study, Taxa, Transect)%>%
  distinct()%>%
  left_join(butterflies_TSI, by = c("Specific.Epithet"))%>%
  left_join(birds_TSI, by = c("Specific.Epithet"))%>%
  mutate(TSI = ifelse(is.na(TSI.x), TSI.y, TSI.x))%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet, TSI)%>%
  summarise(CSI.T = mean(TSI, na.rm = T))
# summary
CSI.T%>%
  group_by(Taxa, Study)%>%
  skimr::skim(CSI.T)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

Variable type: numeric

Taxa Study mean sd
Birds 1998 - 2000 2.59 0.16
Birds 2016 - 2017 2.53 0.09
Butterflies 1998 - 2000 5.34 0.39
Butterflies 2016 - 2017 5.07 0.24 
# t - test

CSI.T%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(CSI.T ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()
Taxa 1998-2001 2016-2017 parameter statistic p.value
Birds 2.588009 2.531490 9.462567 0.8103492 0.4376441
Butterflies 5.338563 5.065191 10.056310 1.5727722 0.1466765 
CSI.T%>%
  ggplot(aes(Taxa, CSI.T, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Community specialisation index \n (Trophic)")+
  scale_fill_brewer(palette = "Greys")+
  scale_y_sqrt()

Niche width of the community did not change over the two studies.

6.2 Habitat Niche

HSI <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  group_by(Study, Taxa, Habitat, Specific.Epithet)%>%
  # calculating abundance of each species in each habitat type
  summarise(n = sum(Abundance, na.rm = T))%>%
  group_by(Study, Taxa)%>%
  mutate(N = sum(n), # total number of indviduals encountered 
         rel.prop = n/N)%>% # relative proportion of species
  group_by(Study, Specific.Epithet)%>%
  summarise(HSI = sd(rel.prop)/mean(rel.prop))%>%
  replace_na(replace = list(HSI = 0))

CSI.H <- abn%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
  group_by(Study, Taxa, Transect)%>%
  distinct()%>%
  left_join(HSI, by = c("Study", "Specific.Epithet"))%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet, HSI)%>%
  summarise(CSI.H = mean(HSI, na.rm = T))
# summary
CSI.H%>%
  group_by(Taxa, Study)%>%
  skimr::skim(CSI.H)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

Variable type: numeric

Taxa Study mean sd
Birds 1998 - 2000 0.54 0.10
Birds 2016 - 2017 0.49 0.02
Butterflies 1998 - 2000 0.61 0.09
Butterflies 2016 - 2017 0.57 0.03 
# t - test

CSI.H%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(CSI.H ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()
Taxa 1998-2001 2016-2017 parameter statistic p.value
Birds 0.5363977 0.4942803 6.719090 1.068760 0.3220632
Butterflies 0.6110527 0.5692785 7.782004 1.210166 0.2616879 
CSI.H%>%
  ggplot(aes(Taxa, CSI.H, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Community specialisation index \n (Habitat)")+
  scale_fill_brewer(palette = "Greys")

Habitat niche of the community did not change over the study periods.

7 Seasonal changes in relative abundance

season <- abn%>%
  mutate(Month = month(dmy(Date), label = T))%>%
  group_by(Study, Month, Taxa)%>%
  summarise(n = sum(Abundance, na.rm=T))%>%
  group_by(Study, Taxa)%>%
  mutate(N = sum(n),
         rel.prop = n/N)

season%>%
  ggplot(aes(Month, rel.prop, col = Taxa, shape = Taxa))+
  geom_point(size = 3)+
  geom_path(aes(group = Taxa), size = 1, linetype = "dashed")+
  scale_color_brewer(palette = "Set1")+
  facet_wrap(~Study, ncol = 1)

8 Site Map

library(raster)
library(sf)

# Raw point data

points <- read.csv("./Data/sites.csv")%>%
  drop_na()%>%
  dplyr::select(-Location)

# Converting to SF points

points_sf <- st_as_sf(points, coords = c('Lon', 'Lat'))

st_crs(points_sf) <- 4326

# Making transect lines

lines <- points_sf%>%
  mutate(Transect = as.factor(Transect))%>%
  group_by(Transect)%>%
  summarise()%>%
  st_cast("LINESTRING")
# getting study extent

study_ext <- st_bbox(points_sf)

ext <- st_as_sfc(study_ext)

# getting osm map

library(osmdata)

# Tamhini WLS

tamhini_raw <- opq(bbox = study_ext, timeout = 100)%>%
  add_osm_feature(key = 'boundary', value = 'protected_area')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
tamhini <- tamhini_raw$osm_multipolygons

st_crs(tamhini) <- 4326

# Roads

tamhini_area <- st_bbox(tamhini)

roads_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'highway')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
roads <- roads_raw$osm_lines

st_crs(roads) <- 4326

# Landuse

landuse_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'natural')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
forest <- landuse_raw$osm_polygons%>%
  filter(natural == "wood")

water <- landuse_raw$osm_multipolygons

water_small <- landuse_raw$osm_polygons%>%
  filter(natural == "water")

st_crs(forest) <- 4326

# waterways

waterway_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'waterway')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
waterways <- waterway_raw$osm_lines

st_crs(waterways) <- 4326

# elevation 

library(elevatr)

tam_elev <- elevatr::get_elev_raster(locations = tamhini, z = 8)

tam_contours <- raster::rasterToContour(tam_elev, levels = seq(0, 1000, 100))
# making final map

library(gridExtra)

fig1 <- grid.arrange(fig1a, fig1b, fig1c, legend, layout_matrix = rbind(c(1, 1, 1),
                                                          c(1, 1, 1),
                                                          c(2, 3, 4)))

ggsave(fig1, filename = "./Figures and Tables/figure1.png", height = 8, width = 8)